home *** CD-ROM | disk | FTP | other *** search
- ;;; -*- Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
- ;;; The data in this file contains enhancments. ;;;;;
- ;;; ;;;;;
- ;;; Copyright (c) 1984,1987 by William Schelter,University of Texas ;;;;;
- ;;; All rights reserved ;;;;;
- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
- ;;; (c) Copyright 1982 Massachusetts Institute of Technology ;;;
- ;;; Original code by CFFK. Modified to interface correctly with TRANSL ;;;
- ;;; and the rest of macsyma by GJC ;;;
- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
-
- (in-package "MAXIMA")
- (macsyma-module rombrg)
- (load-macsyma-macros transm numerm)
-
- (declare-top(special user-timesofar))
-
- ;;; the following code if for historical frame of reference.
- ;;;(defun fmeval3 (x1)
- ;;; (cond ((integerp (setq x1 (meval x1))) (float x1))
- ;;; ((floatp x1) x1)
- ;;; (t (displa x1) (error '|not floating point|))))
- ;;;
- ;;;(defun qeval3 (y1 x1 z)
- ;;; (cond (x1 (fmeval3 (list '($ev) y1 (list '(mequal) x1 z) '$numer)))
- ;;; (t (funcall y1 z))))
-
- (DEFMVAR $ROMBERGIT 11. "the maximum number of iterations" FIXNUM)
- (DEFMVAR $ROMBERGMIN 0. "the minimum number of iterations" FIXNUM)
- (DEFMVAR $ROMBERGTOL 1.0e-4 "the relative tolerance of error" FLONUM)
- (DEFMVAR $ROMBERGABS 0.0 "the absolute tolerance of error" FLONUM)
- (DEFMVAR $ROMBERGIT_USED 0 "the number of iterations actually used." FIXNUM)
-
- (DEFVAR ROMB-PRINT NIL ); " For ^]"
-
-
- (defun $ROMBERG_SUBR (FUNCTION LEFT RIGHT)
- (BIND-TRAMP1$
- F FUNCTION
- (LET ((A (FLOAT LEFT))
- (B (FLOAT RIGHT))
- (X 0.0)
- (TT (make-array $rombergit :element-type 'double-float))
- (RR (make-array $rombergit :element-type 'double-float))
- (USER-TIMESOFAR (cons 'romb-timesofar user-timesofar))
- (ROMB-PRINT NIL))
- (setq X (-$ B A))
- (SETF (AREF$ TT 0)
- (*$ x (+$ (FCALL$ F b) (FCALL$ F a)) 0.5))
- (SETF (AREF$ RR 0.)
- (*$ x (FCALL$ F (*$ (+$ b a) 0.5))))
- (do ((l 1. (f1+ l)) (m 4. (f* m 2.)) (y 0.0) (z 0.0) (cerr 0.0))
- ((= l $rombergit)
- (MERROR "ROMBERG failed to converge"))
- (DECLARE (FLONUM Y Z CERR)
- (FIXNUM L M))
- (setq y (float m) z (//$ x y))
- (SETF (AREF$ TT L) (*$ (+$ (AREF$ tt (f1- l))
- (AREF$ rr (f1- l))) 0.5))
- (SETF (AREF$ RR L) 0.0)
- (do ((i 1. (f+ i 2.)))
- ((> i m))
- (COND (ROMB-PRINT
- (SETQ ROMB-PRINT NIL) ;^] magic.
- (MTELL "Romberg: ~A iterations; last error =~A;~
- calculating F(~A)."
- I
- CERR
- (+$ (*$ z (float i)) a))))
- (SETF (AREF$ RR L) (+$ (FCALL$ F (+$ (*$ z (float i)) a))
- (AREF$ rr l))))
- (SETF (AREF$ RR L) (*$ z (AREF$ rr l) 2.0))
- (setq y 0.0)
- (do ((k l (f1- k))) ((= k 0.))
- (DECLARE (FIXNUM K))
- (setq y (+$ (*$ y 4.0) 3.0))
- (SETF (AREF$ TT (f1- K))
- (+$ (//$ (-$ (AREF$ tt k)
- (AREF$ tt (f1- k))) y)
- (AREF$ tt k)))
- (SETF (AREF$ RR (f1- K))
- (+$ (//$ (-$ (AREF$ rr k)
- (AREF$ rr (f1- k))) y)
- (AREF$ rr k))))
- (setq y (*$ (+$ (AREF$ tt 0.)
- (AREF$ rr 0.)) 0.5))
- ;;; this is the WIN condition test.
- (cond ((and
- (or (not
- (< $rombergabs
- (setq cerr
- (abs (-$ (AREF$ tt 0.)
- (AREF$ rr 0.))))))
- (not (< $rombergtol
- ;; cerr = "calculated error"; used for ^]
- (setq cerr (//$ cerr
- (cond ((= y 0.0) 1.0)
- (t (abs y))))))))
- (> l $rombergmin))
- (SETQ $ROMBERGIT_USED L)
- #+maclisp
- (progn (*rearray tt) (*rearray rr))
- (return y)))))))
-
-
-
- (defun romb-timesofar () (setq romb-print t)) ;^] function.
- ;;; Making the ^] scheme work through this special variable makes
- ;;; it possible to avoid various timing screws and having to have
- ;;; special variables for communication between the interrupt and MP
- ;;; function. On the other hand, it may make it more difficult to
- ;;; have multiple reports (double integrals etc.).
-
-
- ;;; TRANSL SUPPORT.
-
- (DEFPROP $ROMBERG_SUBR $FLOAT FUNCTION-MODE)
-
- (DEFUN ROMBERG-MACRO (FORM TRANSLATEP)
- (SETQ FORM (CDR FORM))
- (COND ((= (LENGTH FORM) 3)
- (COND (TRANSLATEP
- `(($ROMBERG_SUBR) ,@FORM))
- (T
- `((MPROG) ((MLIST) ((MSETQ) $NUMER T) ((MSETQ) $%ENUMER T))
- (($ROMBERG_SUBR) ,@FORM)))))
- ((= (LENGTH FORM) 4)
- (LET (((EXP VAR . BNDS) FORM))
- (COND (TRANSLATEP
- `(($ROMBERG_SUBR)
- ((LAMBDA-I) ((MLIST) ,VAR)
- (($MODEDECLARE) ,VAR $FLOAT)
- ,EXP)
- ,@BNDS))
- (T
- `((MPROG) ((MLIST) ((MSETQ) $NUMER T) ((MSETQ) $%ENUMER T))
- (($ROMBERG_SUBR)
- ((LAMBDA) ((MLIST) ,VAR) ,EXP)
- ,@BNDS))))))
- (T
- (WNA-ERR '$ROMBERG))))
-
- (DEFMSPEC $ROMBERG (FORM)
- (MEVAL (ROMBERG-MACRO FORM NIL)))
-
- (def-translate-property $ROMBERG (FORM)
- (LET (($TR_NUMER T))
- (TRANSLATE (ROMBERG-MACRO FORM T))))
-